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Abstract 

In this work, we discuss some points relevant for stochastic modelling 
of one- and two-phase turbulent flows. In the framework of stochastic 
modelling, also referred to PDF approach, we propose a new Langevin 
model including all viscosity effects and thus that is consistent with vis- 
cous Navier-Stokes equations. In the second part of the work, we show 
how to develop a second-order unconditionally stable numerical scheme 
for the stochastic equations proposed. Accuracy and consistency of the 
numerical scheme is demonstrated analytically. In the last part of the 
work, we study the fluid flow in a channel flow with the proposed viscous 
method. A peculiar approach is chosen: the flow is solved with a Eule- 
rian method and after with the Lagrangian model proposed which uses 
some of the Eulerian quantities. In this way attention is devoted to the 
issue of consistency in hybrid Eulerian/Lagrangian methods. It is shown 
that the coupling is important indeed and that to couple the Lagrangian 
model to an Eulerian one which is not consistent with the same turbu- 
lence physics leads to large errors. This part of the work complements a 
recent article [Chibbaro and Minier International Journal of Multiphase 
flows submitted (arXiv:0912.2045J ]. 

1 Introduction 

Stochastic modeling approaches (also referred to as ProbabiUty Density Func- 
tion (PDF) methods), once developed in statistical mechanics and solid state 
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physics [I], have been shown to be very powerful for the study of turbulent flu- 
ids in presence of complex physics, in particular, for turbulent combustion flows 
[5] and for polydispersed two-phase flows [3J. The main advantage of the PDF 
approach over conventional moment-closure methods is its ability to reproduce 
convection and non-linear source terms without approximations. Moreover, the 
information available by using such methods is the complete PDF (though mod- 
eled) and, thus, all the statistical moments which are related. 

In PDF methods, turbulent closure is achieved through a modeled trans- 
port equation for the joint PDF of some chosen variables, which constitute the 
state vector of the process. In this work (and almost ever) the one-particle La- 
grangian PDF is considered, from which the one-point one-time joint Eulerian 
PDF can be extracted. The resulting modeled one-particle PDF transport equa- 
tion is a Fokker-Planck equation [l], which is a high dimensional scalar equation. 
Thus, traditional numerical techniques such as finite- volume and finite-difference 
methods are possible, in theory, but are not suitable to solve PDF equations in 
practice, since the computational cost increases exponentionally with the num- 
ber of dimensions. On the other hand, Fokker-Planck equations are equivalent 
(in a weak sense) to a set of stochastic differential equations (SDKs). In this 
case, the PDF is represented by an ensemble of Lagrangian stochastic particles 
whose properties are driven by the model SDEs. These stochastic particles can 
be regarded as samples of the PDF and following these particles in time rep- 
resents a dynamical Monte Carlo method, whose computational cost increases 
only linearly with the number of particles. For this reason, the Monte Carlo 
method (or particle stochastic approach) is usually chosen to solve high dimen- 
sional equations and, in particular, PDF equations. 

In turbulent two-phase flows, the SDEs equations of the model contain sev- 
eral mean flelds and they have the general form 

dZ.{t) = Mt, Z, (/(Z))) dt + J2 B.j{t, Z, (/(Z))) dW,{t), (1) 
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where the ( ) operator stands for the mathematical expectation. One way to 
compute mean fields is to extract them directly from the particle properties 
with the help of a mesh, this is the stand-alone method. This method is fully 
consistent but suffers, nevertheless, from some drawbacks mainly due to the 
statistical fluctuations in the particle mean flelds^. On the other hand, hybrid 
methods, where the mean fields are provided by a coupled different numerical 
method, are also possible. These hybrid methods are based on particle-mesh 
techniques [S]: the mean-field equations are solved on a mesh by standard dis- 
cretization techniques whereas the dynamics of the particles are still obtained 
by the time integration of the stochastic differential equations. The main aim 
of such hybrid methods is to improve the efficiency of numerical simulations 
without any important loss of accuracy. Indeed, those methods try to conjugate 
the advantages of moment approach, which provides mean-fields free of statis- 
tical errors and low computational costs, with those of PDF one, which are the 
accurate treatment of some non-linear phenomena and the more detailed level 
of information. 



2 



Recently, the possibility to couple numerical methods of different kind has 
emerged in several scientific domains, for example the kinetic approach with 
the fluid one in DSMC simulations [S]. In fluid mechanics, different strategies 
have been explored to couple PDF methods with other approaches, in partic- 
ular, Moments/PDF [fllSj, PDF/SPH (smooth particle hydrodynamics) [5] and 
LES/PDF [10]. In this work, we are dealing with two-phase flows and we choose 
the framework of hybrid Moments/PDF approach, where, on the one hand, the 
fluid flow is approached by classical Eulerian moment method and, on the other 
hand, the dynamics of the solid particles is directly simulated by a stochastic 
process of the form ([T]). As mentioned above, the mean fields present in the 
stochastic model are provided by the Eulerian solver. The present work com- 
plements a recent article [TT] and thus considers the same configuration treated 
there that is we limit ourselves to a fluid-fluid conflguration, that is, we are con- 
sidering tracer particles with zero diameter and, as a consequence, zero inertia. 
This represents an asymptotic limit case of the general two-phase flow conflg- 
uration, but numerically it preserves the same characteristics concerning the 
exchange of variables between the different methods and, thus, it gives a privi- 
leged position to address issues about the consistency of the method, which is a 
major point for such approaches and which can be investigated with difficulty in 
more complex situations. Indeed, in the Eulerian solver, classical second-order 
turbulent models are generally used, while for the solution of SDEs different 
techniques must be used, as explained later. In this framework, several mean 
flelds can be computed as duplicate flelds in the FV and particle algorithm, 
which raises questions of consistency. From a numerical point of view, we have 
chosen a turbulent channel flow as test-case which represents a classical engi- 
neering situation. 

Three main purposes characterize the present work. The first one is to 
propose a new viscous PDF model which takes into account the viscous terms 
that are present in Navier-Stokes equations. This kind of model may be useful 
when simulations of low- Reynolds number or wall-bounded flows are taken under 
consideration, since in those cases viscous terms are normally not negligible. 
The viscous model proposed is shown to be consistent with exact first and 
second-order moment equations for velocity. The second purpose is to propose 
a general and rigorous framework in which to develop numerical schemes for 
SDEs of the very general form ([T]). In particular, we develop an efficient, stable 
and accurate numerical schemes for the integration of the trajectories of the 
stochastic process. This implies two main difficulties. 

(a) The first one arises from the nature of the stochastic models. SDEs do not 
obey the rules of classical differential calculus and one has to rely on the 
theory of stochastic processes ^2J. In the present paper, Ito's calculus is 
adopted and therefore all SDEs are written in the ltd sense. For stochastic 
processes, several convergence modes are possible only weak convergence 
is under consideration, since the purpose is just to estimate mean quanti- 
ties. A discrete approximation Yt {T stands for a given stopping time) 
converges in the weak sense with order p. if for any polynomial 5, there 
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exists a constant C, function of T, such that 

|(5(ZT)>-(g(YT))KC(r)AiP. (2) 

Due to the mathematical definition of Ito's integral, numerical schemes 
developed for PDEs cannot be applied directly. 

(b) The second difficulty is related to physical constraints. As suggested in 
Minier [T^, the general stochastic model used to simulate general turbu- 
lent flows contains several characteristic time-scales. When some of these 
time-scales become negligible (the system of SDEs is then stiff), various 
sub-systems of stochastic differential equations can be extracted. In other 
words, simplified stochastic models can be obtained from the general one. 
Our second objective is to put forward numerical schemes that can be still 
applicable, and that remain accurate, when the different time-scales go to 
zero |il4u8j. Thus, it should be evident that this corresponds to a practical 
concern, for in the numerical simulation of a complex flow, the time-scales 
may be negligible in some areas of the flows. We want nevertheless the 
general numerical scheme to reproduce the correct physical behavior in 
these areas with the same numerical efficiency. To overlook this point 
can cause severe problems and even flawed results. For instance, in many 
studies, wall-bounded flows are considered where SDEs become stiff near 
the wall. In these cases to use standard numerical schemes like standard 
euler is not stable. The consequence is that very small time-steps have 
to be used to stabilise the numerical scheme. That increases much the 
computational cost, on one hand, and makes results questionable, on the 
other hand. 

The third purpose of this work is to investigate numerically the behavior of the 
method, even in the asymptotic limits, with the main attention devoted to the 
issue of consistency. Indeed, hybrid methods are very attractive for simulations 
of reactive and multiphase flows, but the influence of the exchange of variables 
remains to be completely understood. It is emphasized here that it is possible 
to assess the global consistency of the method, because we use a completely 
consistent and an accurate numerical scheme, which allows us to consider the 
numerical errors as negligible. Moreover, we will deal with the effect of using 
not-consistent Eulerian models. Fianlly, a comparison between viscous and 
high- Reynolds results is carried out. 

2 Viscous Model 

In this section, we propose a new Langevin model which account for the viscous 
terms present in the exact PDF equation. Indeed, starting from the exact 
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equations for a fluid, Navier-Stokes equations 



^ + U-^ = --— + iy^ (3b) 
dt ^ dxj pf dxi dx'j 

we can derive an exact (yet not closed) equation for the one-particle PDF equa- 
tion for the state-vector Y 

-^[(.Af/.IY.y)/], 

this equation may be considered as the starting point for stochastic modeling. In 
high-Reynolds number flows, all viscous terms, but the turbulent dissipation, are 
generally neglected; however, in low-Reynolds number flows or in wall-bounded 
flows, they can become necessary. In order to include those cases also in the PDF 
approach, two different Langevin viscous models have been already proposed 
and tested [121 [H] . Here we propose a third possible model for the state- vector 
(x,U): 

dxi = Ui dt (5) 

dU^ = ^dt + V ^ dt 

p oxi oxj, 

+Aa{Ui - {Ui))dt + Ga{Ui - {Ui))dt 

+ ^Co{e)dW,, (6) 

where the matrix Aij is defined implicetely through the equation: 



I' 



(AR).,, ^ - [^-^j , with R,, = {u^u^y, (7) 

and the matrix dj is the general matrix that models the pressure-fluctuation 
term. In this work, we use the simplifled Langevin model (SLM), that is 

At last, W is a Wiener process [Tj. 

This model has been developed in the same framework of the two others 
cited above, thus, for the general discussion of physical properties of this class 
of models we refer to those papers. In the appendix |X1 we show that the model 
proposed ([S])-® gives the correct equations for the moments of order 1 (mean- 
velocity), and 2 (Reynolds-stress) and, therefore, that the model is consistent 
with the known physics of fluid turbulence. 
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The main specificity of this new model is that the equation ([S]) for the posi- 
tion of the stochastic particles is maintained equal to the exact equation of fluid 
particle-position (in a Lagrangian sense). On the contrary, both other models 
used the artifice of adding a white noise in the modeled equation of particle 
position to represent the effect of diffusion. Whilst this feature is only of formal 
relevance in fluid turbulence, it can be useful in the case of two-phase flows |17) . 
where only this new model can be straightforwardly used. Since the equation 
for particle position is constrained to remain not-modified, an ad-hoc term has 
been added in the modeled equation for particle velocity, such that the diffusion 
in space for mean-velocity and the term constituted by the matrix Aij is nec- 
essary to assure the correct equations for the Reynolds-stress. In general flows, 
the matrix Aij can be computed from the formula ([7]) , that can be rewritten in 
the form 



In present work, we consider the simple case of a channel flow, and further 
simpliflcations can be made considering the symmetry properties of Reynolds- 
stress in X and z directions. Furthermore, we make the hypothesis that the 
Reynolds-stress are nearly constant in the logarithmic zone, which is well verified 
experimentally. This is, in fact, an equilibrium hypothesis and, on his basis, we 
can take just the leading terms in the viscous zone and solve the equation 
analytically. Then, we have 



A • R = e ; hence A = e • R 



(9) 





(10) 



Al2 









^33 



(11) 



with 




Al2 



11 ~ TT 



(12) 



(13) 



(14) 



3 Numerical Scheme 



Given the physical model, two main issues remain: 

(i) Numerical integration scheme 

(ii) Bomidary conditions for particles 



6 



3.1 Stochastic differential system 

The set of SDEs to be integrated can be written in a compact form 
dxi — Ui dt, 

dU^ ^^^U^dt + Ai2U2Si,,dt + C^dt + Y^ a,j dWj {t), 



(15) 



T,, 



where the diffusion matrix is diagonal and is defined 

ay(t,x,(U)) = (e)Co<5y, (16) 
and vector C is the sum of several terms 

C, = M_,4„([/,)-i^ + „v^((/,) (17) 
C, - M_i£^,„v=(a,) (18) 

a . M-i£^ + „v^(a,), (19) 

Ta p dz 
and the time-scales have become 

We recall that this model has a physical meaning only in the case where 
r] ^ dt ^ Ti, with rj the Kolmogorov time-scale. When this condition is not 
satisfied, it is possible to show that, in the continuous sense (time and all 
coefficients are continuous functions which can go to zero) , the system converges 
towards several limit systems [13]: case 1 : when — > with no condition on 
o'ijTi, the velocity is no longer random and the system becomes deterministic. 
The flow is laminar and it can be proven that 

, , 1 dxi = Uidt , , 

system > I , ^ 21) 

^ T.^O \ U, = {U,). ^ ' 

That is the expected behavior in the viscous sub-layer. 

case 2 : When — and at the same time CTyT; est and Ai,2 = est, the 
fluid velocity becomes a fast variable. It can be shown that 

system ^ i'^.jT^'^cst)^ ^ _^ '^(cTijTi) dWj {t). (22) 

i 

We retrieve a pure diffusive behavior, that is the equations of the Brownian 
motion. Actually, this limit is formal in the continuous sense, being impossible 
to fit both conditions in a turbulent bounded flow {a ^ ^/e which is always 
bounded). Nevertheless, in numerical computations this limit can be recovered 
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|14| and, furthermore, it has a theoretical relevance, since those conditions are 
required for the elimination of a variable, treated as fast variable. The same 
reasoning has been applied to the fluid acceleration to obtain the present model 
on velocity Thus, it is important that the numerical scheme is consistent 
even with this formal limit. 

In our model, the diffusion matrix is diagonal, Eq. (fTO)) . Although following 
manipulations are strictly valid for the general case, for the sake of simplicity, 
from now the diagonal expression will be retained. 



3.2 Analytical solutions 

The construction of the numerical schemes is now slightly anticipated. Since 
the numerical methods are derived by freezing the coefficients on the integration 
intervals, the solutions to system (|f 5p . with constant coefficients, are now given. 

The method of the constant variation is used. For the fluid velocity with 
i 7^ 1, if one seeks a solution of the form Ui{t) = H{t) exp{—t/Ti), where H{t) 
is a stochastic process, Ito's calculus gives 

dH{t) = exp(VT,)[C,;(x) dt + cr,(x) dWj{t)l (23) 

that is, by integration on a time interval [^o, t] {At ^ t — to), 

U^{t) = U^{to) exp{^At/T,) + C,(x) T, [1 - exp{~At/T,)] 

+(T,(x)exp(-VT,) f exp{s/T,)dW,{s). 

J to 

We can define 



(24) 



7,(i) = a,(x) exp(-VT,) / exp(s/T,) dW^{s) ; for z ^ 1 (25) 

Jto 

For i — 1 we substitute this expression to C/2 and we proceed in the same 
way to obtain 

Ui{t) = Uiito) exp(~A</ri) + Ci(x) Ti [f - exp(-At/Ti)] 

+ A12 eM-t/Ti) I exp(s/ri) [(73(^0) exp(-As/r2) + C2(x) T2 (f - exp(-As/r2)) + 72(5)] ds 

Jto 

+ (7i(x)exp(-t/ri) / exp(s/Ti)dH/i(s). 

Jto 

(26) 

Performing integration by parts we, finally, obtain 

Ui{t) = ;7i(to)cxp(-A</ri + [Ci(x)Ti+C2(x)T2Ai2r2](f-exp(-At/Ti)) 
+ Ai2(?[t/2(to)-C2(x)T2][exp(-At/r2)-cxp(-At/ri)]+7i(t) (27) 
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having defined 



71 (i) = Ai2exp(-t/ri) / exp(s/Ti)72(s)ds + aiexp(-Vri) / exp{s/Ti)4m} 

Jto Jto 



TIT2 



(29) 



The first integral in the definition of 71 is a double stochastic integral which can 
be simplified through the integration by parts, as explained in appendix IB. II 
The new formula is 



71 W 



exp(-t/T2) f exp(s/T2)dW'2 - exp(-i/ri) f exp(s/Ti)dW^2 



+ CTiCxp(-t/ri) f cxp{s/Ti)dWi 



(30) 



The analytical solutions for the positions are worked out is a similar manner 
starting from Eq. 



dxi = Uidt Xi{t) = Xi{to) + / Ui{s)ds 



(31) 



for i ^ 1 



x^{t) = x,(to) + U^{tQ)T,[l - exp(-At/T,)] + [C,(x)r,]{At - T,[l - cxp(-AVT,)]} + 11^2) 



T,{t) = G, \ &xp{-t'lT,) \ exp(s/T,)dW,(s) 



(33) 



*0 



As usual, the double integral can be transformed 



dW,[s)-eM-t/T^) / eMs/T,)dW,{s) 



to 



; for i = {2, 3} . 

(34) 



Finally, the integration of the first component in position gives 

xi{t) = a;i(to) + f/i(io)Ti[l-exp(-At/Ti)]+Ti[Ci(x)+Ai2C2(x)T2] X 
X {At - Ti[l - cxp(-AVri)]} + Ai2e[U2{to) - C2(x)r2] X 
X [r2(l-cxp(-At/r2))-Ti(l-exp(-AVTi))]+ri(t) (35) 

Ti{t) = Ai2<J2e{{T2-Ti) f dW2{s)-T2eM~t/T2) f exp{s/T2)dW2{s) 

I Jto Jto 

+ T^eM-t/Ti) exp(s/ri)dM^2(s)| 

/ dWi-Tiexp{-t/Ti) [ exp{s/Ti)dWi{s) ; (36) 

Jtn Jto J 



+ CTl 
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also for the position the stochastic integral Fi has been reduced in his basic 
components through integration by parts. To conclude this section, we resume 
all the analytical solutions in table [TJ 

The stochastic integrals, Eqs ([5^ to ([55]) in Table [1] only for the constant 
coefficient case, are Gaussian processes since they are stochastic integrals of 
deterministic functions [T^]. These integrals are quite involved, but. actually, 
they are expressed in function of some simple stochastic integrals which can be 
used as base. The base is formed by the following 7 integrals 



/i,, = exp{~t/T,) exp{s/T,)dW,{s) i = {l,2,3} 



i2,2 



^3,1 



eM~t/Ti) / eMs/Ti)dW2{s) 

J to 

dW^{s) i = {1,2,3}. 

I to 



(37) 
(38) 
(39) 



The original stochastic integrals eqs ([5^ to (155)) are, then, decomposed and 
reformulated as 



i=l 71 = Ai2(y26[ha- ha] + <^ih,i 



(40) 
(41) 



i^l Ti (JiTi[h^i~ h^i] 

i=l Ti = Ai2a2e[{T2-Ti)h^2 



= M 



3,3/ 



with 



M = 



(42) 



T2ha + Ti/2,2] + ^TiTi[/3,i - /i,il43) 
In a matrix form the two set of integrals are related by the formula 



(44) 



CTl 


Ai2a20 





-^12(72^ 












02 




















0-3 











0-1 Ti 


-72^120-2^ 





TiAi2a20 


(TiTi (r2 - 


- Ti)Ai2a20 








— (T2T2 











0'2T2 






















0'3'73 



(45) 

In order to evaluate and to simulate numerically stochastic integrals, it is 
necessary to compute the moments of integrals and in particular the covariance 
matrix C^. Here, we do not proceed to the complicated calculation of the co- 
variance matrix of the integrals T and 7, but we compute the simpler covariance 
matrix of the base integrals ([37|)-([39ll. Once computed this covariance matrix 
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and with the help of the relations given by Eqs. (P^ - (P5|) . it always possible 
to retrieve the covariance of F and 7 and, then, to simulate them numerically. 
The covariance matrix of the base integrals (|571) - p9l) is obtained using the tech- 
niques explained in appendix IB. 21 and the resulting values are shown in table 



Again, we slightly anticipate the numerics and it is explained how the stochas- 
tic integrals can be calculated (simulated). It has just been shown that these 
integrals are Gaussian processes whose means and variances are know (zero 
mean and covariance matrix given by Table [5]) . The vector composed by the 
seven stochastic integrals, eqs ([37]) to ([59)) . is a vector of Gaussian centered ran- 
dom variables which can be computed by resorting to the simulation of a vector 
composed of independent normal Gaussian variables (zero mean and variance 
equal to one). This technique, which requires the Choleski decomposition of the 
covariance matrix, is displayed in Appendix [C] It is worth recalling that once 
computed the stochastic integrals forming our base, the integrals ([5^ - ([55|) can 
be obtained through the matrix relation ([Uj). 

At last, let us check that the expressions of Table [T] and [5] are consistent with 
the limit cases given, in the continuous sense. 

case 1 : when — > with no condition on at, the flow becomes laminar, 
which means that the system becomes deterministic, see Eq. ([21]) . Once again, 
the results given by Table [T] and [2] are in agreement with Eq. ([2T|) . When 
T, 0, eqs ([15]) to ([5T]) become 



which is the analytical solution to system ([^^ when the coefficients are constant. 
case 2 : When — )■ and at the same time cTiTi — > est, we have 



where Qi is A/'(0, 1) vector (a vector composed of independent normal Gaussian 
random variables) . 

A major conclusion can be derived, in order to have a consistent and stable 
numerical scheme for the integration, even within the asymptotic limits, it is 
necessary to retain the exponential terms in the solution and all the stochastic 
integrals, even the indirect ones, must be explicitly calculated. 



m 



U^{t) = ({/,(i)), 
x,{t)^xdt^) + {U,{t))M, 



(46) 



dxi — (Ui) dt + aiTiVKtQi 



(47) 
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Table 1: Analytical solutions to system psp for time- independent coefficients. 



xi(t) = xi{to) + C/i(to)Ti[l - exp(-At/Ti)] + Ti[C7i(x) + Ai2C2(x)T2]x 
X {At - ri[l - exp(-AVri)]} + Ai29[U2{to) - C2(x)T2] x 
X [T2(l - exp(-At/T2)) - ri(l - exp(-AVri))] + ri(t) (48) 

x^{t) = x,{to) + Udto)T^[l - cxp{-At/T,)] + [C,(x)T,]{Ai - T,[l - exp(-AVr,)]} 
+ ri(t) (49) 

Ui(t) = Ui{to) exp(-At/Ti + [Ci(x) Ti + C2(x) r2Ai2T2] (1 - exp(-Ai/ri)) 
+ A,2e[U2{to) ~ C2(x) T2] (exp(-At/T2 - exp(-At/Ti) + j,{t) (50) 

U^it) = C/,(to) exp(-Ai/rO + a(x) [1 - ejcp{-At/T,)] 

+ a.(x)exp(-t/T,) /" exp(s/T,) dW^,(s) (51) 

J to 



The stochastic mtegrals 7i(i), Ti(t), are given by: 



ri(0 = Ai2t72^ 



(Ta-Ti) /" diy2(s)-T2exp(-i/r2) / exp(s/T2)dW^2(s) 

J to J to 

+ TieM-t/Ti)[ exp(s/ri)dM^2(s) 

"'to 

I dWi~Tiexp{-t/Ti) [ exp{s/Ti)dWi{s) 

Jtn Jtn 



(52) 



to "'to 

dW,{s) - exp{-t/T,) I exp{s/T,)dW,{s) 

to J to 



■ for i = {2, 3} 



(53) 



7i(0 = Ai2exp(-VTi) / exp(s/Ti)72(s)ds + aiexp(-t/Ti) / exp{s/Ti)dWi 

J to J to 

(54) 



7,(t) =a,(x)exp(-Vr,;) / exp(s/T,) dW,(s) 

Jto 



(55) 
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Tahlp ').: Dprivatinn of tVip mvarianrp matrix for mnstant mpffidpnts. 



aMW) = y[l-exp(-2Ai/rO] (56) 

{llS)) = ^[l-eM^2^t/T^)] (57) 

{Ilm=^t (58) 

{h,^h,^) = [1 - exp(-At/T,)] (59) 

{h,2h:2) = [1 - exp(-At/ri) exp(-At/r2)] (60) 

(/2.2/3,2) = Ti [1 - exp(-At/Ti)] (61) 



3.3 Constraints of the numerical schemes 

In the particle-mesh method adopted here, the PDEs for the fluid are first 
solved and then the dynamics of the stochastic particles is computed. Thus, 
the scheme has to be explicit. Furthermore, the time step, which should be the 
same for the integration of the PDEs and the SDEs, is imposed by stability 
conditions required by the Eulerian integration operator. This implies that, 
since there is no possibility to control the time step when integrating the SDEs, 
the numerical scheme has to be unconditionally stable. At last, since particle 
localization in a mesh is CPU demanding, the numerical scheme should minimize 
these operations. The first constraint is 

(i) The numerical scheme must be explicit, stable, of order 2 in time and the 
number of calls to particle localization has to be minimum. 

A practical strategy to fulfill the stability condition is to base the scheme on 
the analytical solutions presented in Section [X^ Indeed, the time step appears 
in decaying exponentials of the type exp(— Ai), which brings unconditional sta- 
bility. Therefore, the second constraint is 

(ii) The numerical scheme must be consistent with the analytical solutions of 
the system when the coefficients are constant. 

(iii) The numerical scheme must be consistent with all limit systems. 
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Table 3: Weak first order scheme (Euler scheme): schl 

Numerical integration of the system : 
i^l +1 ^x'l + AU^+B T^'C^ + 

i=l =x'l+A U''^ + B T^[C1 + ^^2^2"?' - 2"] + C A'l^ [U^ - C'^T^] + T^*, 

i^l U-^' = U: exp{- At/TP) + [Trcnoil - exp(-At/7^")] + 7,", 

i=l C/f +1 = f/f exp(-At/Ti") + ri"[Ci" + a;^2C2"T2"]^ + M - c'iri\E + 7f , 



The coefficients A, C, D and £^ are given by : 
A = - exp(-At/Tf ) 
S = Ai - A 

c-r[r2"(i-exp(-Ai/r2"))-ri"(i-exp(-AVri"))] with 0" = ^ 

^ ^ [1 _ exp(-At/Tr)] 

= A^2r (exp(-At/r2") - exp(-At/Ti")) 

(62) 



For a general and comprehensive discussion about the delicate point of 
asymptotic limit, we refer to the paper where this issue is analyzed in the 
more general case of two-phase flows. The construction of the scheme on the 
analytical solutions should also ensure a sound physical treatment of the multi- 
scale character of the problem. Indeed, it has been demonstrated that from the 
analytical solutions given in Section 13.21 the limit cases can be retrieved. 

3.4 Weak first-order scheme 

The derivation of the weak first order scheme is rather straightforward since 
the analytical solutions of system ([T5|) with constant coefficients have been eval- 
uated. Indeed, the Euler scheme (which is a weak scheme of order 1 [18]) is 
simply obtained by freezing the coefficients at the beginning of the time inter- 
vals At = [t„, t„+i]. Let Z" and 2'"'^^ be the approximated values of Zi at time 
tn and tn-f^i- The Euler scheme is then simply written by using the results of 
Tables [1] and [2] as shown in Table |3l This scheme fulfills all criteria listed in 
Section [3.31 except, of course, the order of convergence in time. 

As far as the consistency with all limit systems is concerned, some precisions 
must be given. Here, it has to be understood that the limit systems are presented 
in the discrete sense. The observation time scale dt has now become the time 
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step At. The physical time scale Ti does not go to zero, as in the continuous 
sense, but its value, depending on the history of the particles, can be smaller or 
greater than At. 

In the limit case 1 (PT|) we have, with 1 ^ At/Ti, 

r,{t) > a, T.VAtg,,, w (63) 

j,{t) > (64) 

and, thus, for the fluid variables 

Ul'+' - m (65) 
<+i = x^ + {UnAt. (66) 

In this case these relations are identical to those analytic. In the second limit 
case it is different; for the present scheme, Eq. (|5ip. with 1 <C At/Ti and 
the constraints on cr, gives 



ur+' = {ur) + ^ (67) 

From the continuous point of view, the fluid velocity Ui will behave as a fast 
variable and a near white-noise term. On the contrary, in the numerical simula- 
tion, the velocity Ui does not, of course, behave strictly-speaking as white-noise 
function (its variance is not infinite!). This result is physically sound. Indeed, 
when a fiuctuating physical process (random variable) is observed at time steps 
which are greater that its memory, the expected behavior is Gaussian. For 
particle position the scheme gives 

x^' =x7 + m)At + ^^ofT^t (68) 

as expected. 



3.5 Weak second order scheme 
3.5.1 Property of the system 

The diffusion matrix of system (jlSp has a singular property, that has a crucial 
importance here [2]. From Eq. (fTC|) . it can be noticed that o-^ depends only 
on time, position and the mean value of the relative velocity. Therefore, the 
only variable of the state vector whose is a function is Xp. Therefore, the 
diffusion matrix has the following singular property 

EE-'^.fe-O' (69) 
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3.5.2 General method 



Let us consider, for a moment, the following stochastic differential equation 

dX,{t) = Ai{X{t)) dt + Y, Bij(X{t)) dWj{t) 
j 

where A is the drift vector and B is the diffusion matrix. If B verifies the 
property (|69p. it can be shown that, for example by stochastic Taylor expansions 
[TB] , a prediction-correction scheme of the type 

3 

is a weak second order scheme (^2^^^^ = ^(taX^+i) andtaS,"^^ = Sy (t2X"+i)) 
This result is true, one again, only when the property is verified. If not, 
other terms are needed in order to obtain a weak second order scheme, see for ex- 
ample Talay [19j. As a consequence, it this scheme is used for a set of stochastic 
equations which do not verify the property (|69p. problems of consistency arise 

m- 

The idea is now to built from Eq. ()70p a weak second order scheme for 
system The first-order Euler scheme is going to be used as a prediction 

step. The remaining task consists in finding a suitable correction step which 
enforces the three conditions listed in Section [3731 



3.5.3 Implementation of the scheme 

The main idea for the implementation of the correction step is to start from 
the analytical solutions to system when the coefficients are constant and 
apply the same idea as in Eq. (j70p considering that the acceleration terms vary 
linearly with time. 

Before the implementation of the numerical scheme, let us specify the nota- 
tion which is in use. As in Section [3741 Z" and Z""''^ stand for the approximated 
values of Zi at time t„ and respectively. Therefore, Z'^'^^ is the corrected 

value whereas the predicted value is denoted <2^"^^, *-e- the value of Zi pre- 
dicted by the Euler scheme. The predicted velocities are t2U"^^ and the values 
of the fields taken at {tn+i, x"+-'^) are denoted, for example, {P'"^-^). The pre- 
dicted time scales are referred to as t2T"^^ and the predicted diffusion matrix 
is designated by t2cr"^^. 

As far as the computation of the fields 'attached' to the particles are con- 
cerned, it is worth emphasizing that none of them are computed at (t„+i, x"'+-'^), 
because the scheme would become implicit. All fields which characterize the 
fluid, such as the mean pressure field, are taken at x"+^), but fields such 

as the expected value of the fluid velocity are computed from the predicted 
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velocities. For example, one has 



-(c/f+i) + n,((u"+i),(p"+i)) 



The analytical solution to system p5)) when the coefficients are constant is, 
for the fluid velocity seen, by applying the rules of Ito's calculus 

i^l U,{t) = C/,(to) cxp(-Ai/T^O + / a{s) exp[(s - t)/T,] ds + j,{t). (71) 



As anticipated we suppose that the acceleration Ci(s,Xp) varies linearly on the 
integration interval [ioj^]j that is 

a;(s,x(s)) = aitoM-to)) + ^[Q(<,x(0) - C,(to,x(io))](s - to). (72) 
By direct integration, one can write 

U^{t) = U,{to) exp(-AVr,)+[r,C,(to,x(io))] A(At,T,)+[T,C,(t,x(t))]B(At,r,)+7^(i), 
where the functions A{At,X) and B{At,X) are given by 

"1 - exp(-AV^)" 



A{At,X) = -cxp(-At/X) 



At/X 



B{At,X) = 1 - 



1 - exp(-Ai/^) 



At/X 



(73) 



By direct application of the ideas presented in the last section, it is proposed to 
write 



C/f+i =i C/r cxp(-AV7^") + i C/r cxp{-At/t2T[^+') 

+A(At, rn [Tl'C^] + B{At, hTl'+^) [t2Tl'+H2C^+^] + t2Y;+\ 



(74) 



The meaning of the corrected stochastic integrals will be given later. For the 
first components the exact solution is 



=1 Ui{t)^Ui{t 



o) exp(-AVri) + / Ci(s) exp[(s - t)/Ti] ds 

J to 



+ Ai2C/2(s)exp[(s-t)/Ti]ds+ / exp[{s~t)/Ti]c7i{s)dWi{s) . 

J to J to 

(75) 

Here, we introduce the exact expression for the variable U2 

U2{t) = U2{to) cxp(-At/r2) + / a(s) exp[(s - t)/T2] ds + 72(t). (76) 

Jto 
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and, then, we let the functions Ci and a varying hnearly. The subsequent 
integration, quite involved, permit us to propose the following expression 

(77) 

where the functions Ac{X,Y) and Bc{X,Y) are given by 

A, {X, Y)^- cxp(- At/X) + [1 - exp(- At/X)] - (^1 + C,{X, Y), 

B, {X, Y) = l- ^^[1 - exp(-At/X)] + 1^ C,{X, Y). 
with 

C, {X,Y) = [exp(-At/r) - exp(~A</X)], 

The stochastic integrals t27i"^^ and t2ji"~^^ remains to be defined. Their 
expression are computed using the same hypothesis, i.e. we suppose that 

(Ti(s,x(s)) = cri(to,x(<o)) + ■^[(Ji{s,x{s)) ~ (Ti{to,x{ta))]{s ~ to). (77) 

After the integration by parts has been carried out, all the stochastic integrals 
are manipulated in order to express them as functions of the base integrals 
([37|) - (j39l) . as done in section |321 Then, final expressions are obtained 

7,"+^ = KA^i2At,t2T['+')+al'+^B^{2At,t2Tl'+')]Ii4t2T['+') i={2,3}(78) 

7r+^ = KA^(2At,<2Ti"+')+<+iB^(2At,t2Tr+^)]/i,i(<27;"+') i=l 

+ <2A"+42r+iKA^(2At,t2r2"+i) + a'^+^B^{2MMT^^^)]h,2{t2T^^^) 

~ <2A"+42r+Vjv4^(2Ai,t2Ti"+i) + a'^+^B^{2At,t2T^+^)]l2,2it2T^t^) 

where the functions A^{At^X) and B^{A.t,X) are given by 

1 



A^{At,X) = 
B^{At,X) = 



/ „A 1 - exp(-2At/X) 

-cxp(-2At/X) + - ' ' 



_^ 1 - cxp(-2At/X) 



2At/X 



2At/X 
1 



1 - exp(-2At/X) 



1 -exp(-2A</X)' 

(80) 



It is essential that the stochastic integrals are simulated by the same Af{0, 1) 
random variable used in the simulation of the stochastic integrals I^j computed 
in the Euler scheme. In appendix l3.6l it is shown by stochastic Taylor expansion 
|18| that the present scheme is a weak order scheme of order 2 in time for 
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system ([T5|) . It is worth emphasizing that no correction is done on position, 
X, since the prediction is aheady of order 2. This property is in line with the 
constraint stated in Section 1X51 that is, the numerical scheme should minimize 
the procedures where particles must be located in a mesh (which is done every 
time the particles are moved, i.e. when a new value of x is computed). The 
complete scheme is summarized in Table 2] 



3.5.4 Limit cases 

When the flow becomes laminar, that is when — ^ with no condition on the 
product (TiTi, one has the following limits: A{At,Ti) — >■ 0, B{At,Ti) — >■ 1 and 
7i(i) 0, which gives for the fluid velocity. 

It can rapidly be shown, by regular Taylor expansion, that this scheme, together 
with the prediction step (Euler scheme schl) is a second order scheme for system 
([2T|) . In the other limit, the fluid velocity become a fast variable (limit case (ii)), 
that is when 1 <^ A.t/Ti and ai Ti est, one can write 

ur' = iur') + (81) 

which is in line with the previous results. In limit case (ii), the numerical 
scheme for the position of the particles is equivalent to the Euler scheme written 
previously and is of first order in time. 



3.6 Consistency and order of the scheme 

In this section, we show by stochastic Taylor expansion that the numerical 
scheme proposed is weak second-order. For a general set of SDEs 

dyi = Di{y)yi + aij{y)dWj 

the general condition to be fulfilled is that [TB] 

y^+' - y^ + D^M + a^^AW, 

1 ^^D?\ 1 rda, 



cT^.AtAW, + hr^ D'^^tAW, (82) 



2\dykJ ' 2\dy, 

+ i;\7^]<'knyir,Atm (84) 



4 \ dykdyi 



da- 



dVk 



WidWj (85) 
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In present case we have 



with ai — 1/Ti, thus 
dVk 



( ^ 


1 




\ 









1 













1 


X3 





-ai 







Ui 












U2 


V 








-as/ 





dykdyi 



— PijSjk H 



dyk 



dykdyi 



dyi 



/aft 



V 9yk 



In the equations given in tables |3]|4l we proceed to Taylor developments of 
all deterministic functions. Moreover, those quantities which are marked by a 
tilde in the correction step are explicitly written and, then, developed. In this 



way, all functions are computed at time n, and Sy = y 



n+l 



y" 



thus, for the 



sake of simplicity, the upper-script is canceled from now on. We obtain for the 
prediction step 



5t2Xi 

5t2Ui 



{-aiyi)At + aiiAWi 

i(aiCTii)AtAVKi + '^{Ai2(J22)AtAW2 
1 



1. 



a2)Ai2At^ 



(86) 



5t2U, = {-alyl+■i)At + (T,iAW^ 

+ i(a,')At' - i(a,a,,)AtAW, with i = {2, 3} . (87) 

Then, for the correction step we use the properties of matrix <j which is diagonal 
(cij = cTiSij) and respect the property (l69l) . For the sake of simplicity we use 
the tensor notation for partial derivatives. 
The stochastic integrals t2"/^~^^ give 

t27r' = (OAt) 

+ + ^{dk<yr)akAW,AWk 0) (88) 

+ ^{dk<y,)PkiyiAW,At + ^{d,l3,,)AWjAtiOAt^/^) (89) 

+ ^idkPrj)a^akAWkAW,At (OAt^) (90) 

+ ^{dk<Jj)l3,jakSk,At^ (=0) (91) 
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finally, we globally obtain 



yf+i = y'l + D,At + a,AW, 

+ ]-AtAW^ [dk(T{)hiyi + (dkMyi'^k + (92) 



2 



Pl^yj + [dkMPkniyiym. + {dkPtkWk + ■7^{dkil3ij)<^k<yiyj^^) 



In the present case, the last expression coincides exactly with the equation re- 
quired by Taylor expansion in the general case and, therefore, it is demonstrated 
that the numerical scheme proposed is of order two in the weak sense. 

In summary, a weak-second order scheme for system (jl5p has been derived. 
This scheme satisfies all conditions listed in Section [231 except for the second 
order convergence condition in limit cases (ii). This imperfection is bearable 
since this case has mainly a formal importance and it does not occur very 
often in practice, moreover, it is inherent to the spirit of the scheme, that is a 
single step to compute position x (in order to minimize the number of particle 
localizations in the algorithm). 

Finally, we want to recall which are the principal hypothesis and proper- 
ties used and what steps have been followed to achieve this weak second-order 
numerical scheme. 



a) The analytical solution to the system p5)) is calculate, taking constant 
the coefficient. It is worth emphasizing that it always possible (almost 
formally) to derive analytical solutions for linear SDEs with constant co- 
efficients [T]. 

b) The complete covariance matrix of all stochastic integrals is calculated. 
In the present case, the stochastic integrals have been preliminary decom- 
posed on a base of simple stochastic integrals and, then, the covariance 
matrix of these last ones is computed. The original integrals appearing in 
the equations can be easily computed through a matrix equation. Gen- 
erally speaking, this operation is not necessary, but sometime it can be a 
powerful method to simplify involved calculations. 

c) The analytical solutions are discretised to obtain an Euler weak first-order 
scheme. Since this scheme is derived directly from the analytical solutions, 
the consistency of the scheme, even in asymptotic limits, is naturally as- 
sured. Moreover, the presence of exponential terms does guarantee the 
global unconditionally stability, though the scheme is explicit. 

d) Thanks to the particular structure of the diffusion matrix a, Eq. (|69p . it 
is possible to chose a predictor-corrector algorithm, with correction only 
on fluid velocity, in order to develop a weak second-order scheme. The 
first-order scheme is taken as prediction step. 

e) The correction is developed starting from analytical solutions of fluid ve- 
locity and letting all coefficients to vary linearly in time-step (hypothesis 
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of constant gradient). The basic idea is that the variation of higher or- 
der ( hnear gradient etc) can contribute with terms at least of order 3 in 
time-step, as all coefficients appear inside integrals. 

It is worth underlying, in conclusion at this section, that the numerical scheme 
proposed remains of weak second-order even in the high-Reynolds number flows, 
that is for = 0, when the model becomes the well known SLM model j21) . 
Furthermore, it ever fills all the constraints asked in section 15751 

4 Boundary conditions 

In the present Lagrangian approach the no-slip and impermeability conditions 
are satisfied by imposing the boundary conditions on stochastic particles. A 
stochastic particle is assumed to have crossed the wall boundary (placed at 
2/ = if at the end of time-step Y(t + At) = yout < 0. For those reflected 
particles we impose the conditions 

y^n = IVoutl ; and U(t + At) = . (94) 

Symmetry conditions are imposed at the other boundary placed a the center of 
the channel y — h/2. Along the x direction, periodic conditions are imposed. 

The problem of boundary conditions for walls in stochastic method remains 
a fascinating and challenging open issue. 

5 Numerical Results 

This section present some numerical results which complement numerical sim- 
ulations discussed elsewhere 

A channel-flow is solved both with an Eulerian method and, then, with 
the PDF one. Being a particle-mesh method, particles are moving in a mesh 
where, at every cell, the mean fields describing the fluid are known. Generally 
speaking, the statistics extracted from the variables attached to the particles, 
which are needed to compute the coefficients of system ([S])-® are not calculated 
for each particle (this would cost too much CPU time) but are evaluated at 
each cell center following a given numerical scheme (averaging operator) . These 
moments can then be evaluated for each particle by projection. This the general 
principle of particle-mesh methods: exchange data between particle and mesh 
points. Once again, the expected value for functionals of the state vector are 
not computed directly for each particle but are evaluated at discrete points in 
space and then calculated for each particle by interpolation. The averaging and 
projecting operators are computed by near grid point (NGP) method [F . In our 
particular case, the same physical case is considered by both methods, thus, first 
the flow is computed by the moment method, then, the computed mean flelds 
are used, frozen in time, in the simulation of stochastic particles. Therefore, 
there is not two-way coupling between the two methods and issues concerning 
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averaging and projection of variables are not much influent, provided that an 
interpolation enough refined is used. It is indeed important to note that the 
interpolation and backward estimation schemes introduce some errors [2^ I23j . 
Nevertheless, for the present simulations, the very fine mesh used in the DNS 
(128 points along y) is also used here, which should assure a good interpolation 
of the mean fields to the Lagrangian stochastic particles and make these errors 
negligible in this case. Furthermore, it is worth saying that results obtained in 
the deterministic context do not generally apply to the stochastic one [S]. 

5.1 Consistency 




Figure 1: (a) Mean velocity Results in Wacawczyk et al. [16J (War04), in the 
present Hybrid configuration and in the DNS. The results are in non-dimensional 
units, (b) Reynolds stress results in War04 and in present Hybrid configuration, 
the results are in non-dimensional units. 

The computations have been performed for the case of fully developed tur- 
bulent channel flow at Re = Urh/v = 395. The channel flow DNS results of 
Moser et al. |24] have been taken for comparison as reference data. For the 
simulations, the mesh used in DNS (128 points along y) is also used here, which 
assures a good interpolation of the mean fields for the Lagrangian stochastic 
particles. Quantities designed with the upper-script + arc non-dimensionalised 
with wall parameters. Given that the PDF used is consistent with the Eulerian 
Rotta Reynolds-stress model (see appendix |A|) . this model should be used in 
the Eulerian part of the Hybrid method in order to check the consistency of 
the Lagrangian part. To ensure the best possible consistency, we use the mean 
fields obtained with an equivalent Lagrangian model [16] in a stand-alone con- 
figuration. In practice, the Eulerian mean-fields computed from the Lagrangian 
stand-alone simulations carried on by Wacawczyk et al. [163 (designed as War04) 
are employed, here, as Eulerian counterpart. Indeed, the equations for the mo- 
ments of first and the second order are exactly the same for both models and, 
thus, a perfect consistency is assured. Nevertheless, in that stand-alone calcula- 
tion turbulent dissipation e is not computed directly, therefore, for this variable 
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DNS values are used and, thus, this remains a possible source of inconsistency. 

In figures [1] three different profiles for the mean-velocity are shown: the 
profiles obtained in the stand-alone configuration in the paper of Wacawczyk 
et al. with an analogous viscous model, the profiles obtained in the hybrid 
configuration and the DNS profile, for physical reference. In figures [U the 
profiles of second-order moments, that is the Reynolds stress tensor, are shown, 
for the present and for the War04 calculations. Generally speaking, since both 
methods face the same test-case, if they were completely consistent, the results 
would be equal both for first and second-order moments. The profiles deriving 
from the present PDF method are in good agreement with those computed in the 
stand-alone configuration. Some small differences are normal for the different 
numerical framework and since a different profile is used for e . 

Yet, in figures [U it is possible to observe that in the zone oi y ^ 0.05 — 
0.1{y^ ~ 10 — 50) the first components of Reynolds stress Ru attain a nearly 
constant level. This merits a particular attention. This behavior can be traced 
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Figure 2: Comparison between dissipation and production. 



back to the value of turbulent dissipation e used in the present calculations. In 
figure [21 the dissipation and production profiles are shown as well as their ratio. 
It is clear that in the zone of interest a relation of equilibrium exists, the value 
of the ratio is near to 1 and, then, it becomes of about 0.8 till y+ ~ 100. In 
this range, the Reynolds-stress equations lead to the relations 



= k 



Co + 2V/e 



k = 



1 + 3Co/2 



1 + 3Co/2 ' 
and, then, using the hypothesis V /e=l 



= U^kr: 



, Cn 



(95) 



(96) 



therefore, it is correct to obtain constant value of the first Reynolds stress com- 
ponents in the range (y+ ^ 10 — 50), in the present configuration. To complete 



24 




Figure 3: (a) Reynolds stress, final comparison, (b) Present results in the hybrid 
configurationj and DNS. 



our numerical validation about the consistency of the method, a slight ad-hoc 
modification of e has been worked out. In practice, the relation between pro- 
duction and dissipation has been imposed to be 0.9 in the region y+ ^ 10 — 50, 
while in the rest of the domain the DNS values are maintained. In figure [31 the 
results are shown. The two profiles overlap. That emphasizes the consistency 
of the numerical method used, since the features of the Lagrangian model are 
perfectly reproduced. We recall that the ingredients necessary to reach this 
objective are: 

(i) Consistent physical model. 

(ii) Consistent numerical scheme 

(iii) Accurate global numerical method , concerning also the exchange of in- 
formation from Eulerian solver to Lagrangian one. 

Given that the method is found to be consistent, present calculations are also 
compared with DNS profiles in figure [1] and [3l From a physical point of view, 
the results are coherent with model features. Indeed, it is well known that the 
SLM model, like the Rotta RSM model, underestimates greatly the value of 
the components Rn and is rather isotropic in the other diagonal components 

imiis]. 

5.2 Asymptotic limits 

A study of the numerical method in some asymptotic cases is necessary to con- 
sider the numerical validation achieved. Two cases are chosen, (a) the laminar 
limit case, essential in wall-bounded flows for the presence of the viscous sub- 
layer, and (b) a academic case where the matrix A goes to infinity, thus the 
time-scales Ti and T3 are taken to be zero, and T2 = T^. 

From a theoretical point of view, in these limit cases the expected behav- 
ior can be anticipated. Starting from Reynolds-stress equations, by using all 
hypothesis valid for the channel flow and considering the stationary case, the 
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Figure 4: (a)Laminar limit case, (b) Theoretical asymptotic limit case. 



Reynolds-stress equations can be simplified. The equations for first, second 
diagonal components and for the shear stress are 

^ - 2G22{v^) + withG22 = -l/ri, ^ («2)«Zla2 ^(97) 

oy 2 

= -2M^ + 2Gii(7/2)+2Ai2M, withGii --l/r{98) 
oy ay 

= -2(«2)^ + 2GnM+2G22M+2Ai2(z;') (99) 
oy oy 

thus, as expected, in the laminar case (a) we have to find 

T, ^ ^ = . (100) 

In the other case(b) we have 

An = A33 = A12 ^ c» Ti = --^ = Ts = --^^ > (101) 

All ^33 

hence, from the equations ([Mjl - dM]) we deduce 



(u^) = —iuv) , (uv) = — (v^) = — — cr? 



(102) 



The numerical results in figure |4] show that the asymptotic limits are per- 
fectly recovered without changing time-step. 
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5.3 Hybrid-method error 

In the last section we have studied and verified the consinstency of our hybrid 
method. In particular, we have checked the behaviour of the method in a 
theoretically consistent configuration and in the asymtotic limits. Since the 
consistency of the method has been successfully proofed, we can now study 
the global error which is eventually introduced by using an hybrid method not 
completely consistent. In order to better explain this very important point, let 
us start by considering the two phases separately. 

The hybrid method for fluid-particle system exploites two different numeri- 
cal (and theoretical) approaches. For fluid phase, a classical moment approach 
is used. This kind of algorythm is characterised by the deterministic discreti- 
sation error that is related to the numerical integration of equations in time 
and space. For particle phase, a Monte Carlo method is used to simulate the 
PDF of the stochastic process which describe particle dynamics. This method is 
affected both by deterministic discretization errors and, further, by the Monte 
Carlo error, which is related to the finite number of stochastic particles used in 
numerical simulations. In the stand-alone approach to stochastic equations for 
turbulent flows, the mean quantities present in models of form ([T]) are computed 
starting from the same stochastic particles. This procedure causes a typical Bias 
error which is found to be dominating in stand-alone methods [3]. In hybrid 
methods, this error is avoided by using a moment approach to compute the 
mean quantities that are necessary. In fact, this is the very strenghtness of 
the method. Nevertheless, even in hybrid configurations the introduction of 
variables computed externally may be source of errors. In particular, it is not 
obvious a priori what happens when the mean quantities are computed by a 
method which is not completely consistent with the PDF one. The error intro- 
duced by this inconsistency is intrisically inherent to the kind of hybrid method 
used. Moreover, from a practical point of view, the best way to investigate 
this effect is in a hybrid fluid-fluid configurations, where the influence of mean 
variables can be easily traced back. 

To isolate the hybrid error, all the other errors have been made negligible in 
the following simulations. In both methods, a time-step of 10"* and a spatial- 
step of 10^"* have been used, which can be considered as zero for our purposes. 
5 * 10"* particles have been employed, which can be considered infinity. To 
the sake of clarity, the configuration discussed in last section, where stand- 
alone mean variables were used, it will be called standard configuration in the 
following. 

First, we consider the following configuration: DNS results are used to pro- 
vided the PDF solver with all necessary mean quantities. In figure [SJ we present 
the mean velocity computed by the PDF method in this configuration. The 
DNS mean velocity is also shown as well as the mean velocity computed by 
PDF method in the consistent PDF-PDF configuration. The mean velocity so 
obtained is perfectly in agreement with DNS one and, thus, is strongly different 
from that obtained in the standard configuration. In fact, the exact profile is 
now recovered. In figure [SId, the Reynolds stress are shown for the same con- 
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Figure 5: (a) DNS-PDF configuration: mean velocity, (b) DNS-PDF configu- 
ration: Reynolds stress. 



figuration. The profilos arc dramatically changed in comparison with standard 
results. The {uv) and, as a consequence as we have seen, the (m^) components 
are strongly overpredicted now. On the contrary, the other diagonal components 
are not changed. 

This behaviour can be explained, at least approxiamatively. We have seen 
in the last section, analysing the behaviour in the asymptotic limits, that the 
diagonal (w^) and {w'^) components are essentially dependent on Lagrangian 
time-scale and on diffusion coefHcient, which have not changed and, therefore, 
they remain the same. The cross shear (uv) depends on the gradient of the mean 
velocity, other than on turbulent kinetic energy and on Lagrangian time-scale, 
thus it is much sensisitive to changes in mean velocity and, in fact, the effect 
of this is enourmous. Indeed, the present mean velocity profile atteins a higher 
maximum than in the standard configuration and, thus, it is much steeper. The 
shape of (u^) is a direct consequence of this behaviour. 
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Figure 6: (a) V2F-PDF configuration: Mean velocity, (b) V2F-PDF configura- 
tion: Reynolds stress. 
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Now, we change again our configuration and we couple the PDF solver to a 
low-Reynolds number RANS model, known for his good performance in bound- 
ary flows, the v2f model [35] • In the present calculations, we use a refinement 
version of the model [55]. In figure [SI we show the mean velocity profile togheter 
with DNS and standard PDF results. As previously, the mean velocity given 
by our PDF method is quite near to the Eulerian mean provided to the PDF 
solver, in this case computed by v2f solver. Furthermore, this result is also 
in good agreement with DNS result. In figure [51 Reynolds stress profiles are 
shown. It can be seen that results are quite similar to those obtained in the 
DNS-PDF configuration. This time, the < uv > and the < > components 
are less overpredicted but they show qualitatively the same behaviour. 

This series of results leads us to some conclusions. 

i) In all three configurations tried, the mean velocity computed by stochas- 
tic particles collapses on the value of the Eulerian mean velocity that is 
provided by the Eulerian solver. This is in line with the physics of the 
model, which is basically base on a return-to-equilibrium idea [27.. 

ii) In all cases, but the standard one, there is a difference between Eulerian 
and Lagrangian results at the level of second-order statistical moments 
(Reynolds stress) . This is a conseguence of coupling in a Hybrid approach 
two methods which are not consistent, this reasoning is valid both for 
PDF-DNS and for PDF-V2F configurations. Thus, it is possible to iden- 
tify the global error of hybrid Eulerian/Lagrangian method by comparing 
the actual Reynolds stress profiles with those computed in the complete 
consistent configuration [3] 

iii) In the last two configurations tried here, the results obtained are quite 
similar in practice, even though DNS and V2F approaches are quite dif- 
ferent from a theoretical point of view. This is reasonable. In fact, the 
profiles provided by the Eulerian solver to the Lagrangian one are mean 
velocity, mean pressure, turbulent energy and turbulent dissipation. For 
these variables, the V2F approach gives results in very good agreement 
with DNS ones. Therefore, from the point of view of the hybrid method 
DNS and V2F approaches are very near. Nevertheless, the results which 
have been obtained using V2F are nearer to those consistent than those 
obtained using DNS, figure [H Therefore, V2F model is found to be more 
consistent with present Lagrangian model than DNS, as expected. 

The large error in the off-diagonal component of the Reynolds stress (uv) 
might cause some doubts about the results shown in this section. It is possible 
to verify by a simple numerical trick that these results could be expected. 

Starting from the mean profiles obtained by the V2F approach, we can com- 
pute the {uv) component of Reynolds stress by using the k — e formula 

{uv)^C,-^. (103) 
e dy 
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Figure 7: Reynolds stress from V2F results through fc — e formula. 



The resulting profile is shown in figure [7] and it is compared with that obtained 
in the hybrid PDF-V2F configuration for the same variable. The profiles are 
qualitatively similar. This shows that the PDF method used is approximatetly 
consistent with the k — e, even though it is rigorously consistent with the Rotta 
Reynolds stress model. Furthermore, it demostrates that the difference between 
the V2f model and the k — e model introduces an error which generates the 
behavior encountered for (uv). 

5.4 Standard model Vs Viscous model 

In this work we have proposed and tested a new Langevin viscous model for 
turbulent fluid flows. Other propositions have been recently made, with the 
aim of improving the present description of near- wall layer in the framework of 
PDF approach. In this section, we compare the viscous model with the standard 
high-Reynolds number one, in the hybrid simulation of turbulent channel flow. 
The high- Reynolds number model is generally used with wall-function boundary 
conditions, without integrating up to the wall. In the following computations, 
we proceed to the up-to-the-wall integration and we use the same boundary 
conditions imosed in the viscous case|4l 

First, we carry out the simulation of the channel flow in the standard con- 
sistent configuration (section l5.ip . In figure [8l the Reynolds stress obtained 
in both cases are shown. The differences are not important, in practice, the 
two models give the same results. In figure [51 the Reynolds stress obtained in 
the hybrid DNS/PDF configuration are shown. Also in this inconsistent con- 
figuration, the results obtained with standard and viscous model do not show 
significant differences. 

In figures 1^1 the probability density functions of the normal velocity are 
shown for two locations: y+ = 5 and ?/+ — 10. The PDFs are computed with 
both standard and viscous PDF methods and, for reference, also DNS results 
are shown. Once again, the new viscous terms bring very small changes to PDF 
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Figure 8: (a) Reynolds stress for standard and viscous model, in the consis- 
tent configuration, (b) Reynolds stress for standard and viscous model, in the 
DNS/PDF configuration. 




Figure 9: (a) Probability density functions computed by viscous, standard PDF 
and DNS approach at y"'' = 5. (b) Probability density functions computed by 
viscous, standard PDF and DNS approach at y+ — 10. 



profiles. From a general point of view, it is possible to observe that in this 
configuration with the SLM model the PDFs remain quasi-Gaussian and the 
distorsion due to walls is scarcely felt, mainly in the transition zone at = 10. 
In that sense, given the good results shown by PDF methods equipped by elliptic 
relaxion model [28l [29] , it is possible to conclude that in order to improve the 
quality of the physical picture the essential is given by non-local representation 
of pressure fluctuating term. 

In conclusions, the viscous correction in near-wall models are important 
from a theoretical point of view, for they lead to the exact expression for first 
and second moment equations. However, from a practical point of view, they 
are negligible in hybrid configurations. Indeed, the numerical results are not 
influenced in a noticeable way by the viscous corrections. Other viscous models 
have been proposed [28[ [T5] . but this is the first time that such comparison 
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has been made. Moreover, the other models have been studied in stand-alone 
configurations, where the integration up to the wall is more delicate. Thus, 
it is not possible to generalize this conclusion to all PDF methods, in stand- 
alone methods viscous terms may be necessary in order to compute directly the 
near-wall zone. 

6 Conclusions 

First, a new physically-consistent PDF model is proposed, then, a numerical 
scheme for the stochastic differential equations which appear is detailed. Fi- 
nally, a numerical validation of the global method is carried out with the main 
attention pointed on the issue of the consistency between the two approaches. 

The model proposed is in the form of a set of stochastic differential equations 
©-(in]). Since the model belongs to the same category of those studied previously 
fl5i il6j , a physical validation of the model, to be done necessarily with a stand- 
alone code, should not show particular new insights. Therefore, we suppose that 
the performances of the present model are equivalent to those already validated 
[TB] and we have chosen to put the emphasis on the numerics. With regard 
to thos previous models, the present one is consistent with the same first and 
second order moment equations and thus represents just an alternative. Its 
main feature is that takes into account of the viscosity without changing the 
x-equation, as done in previous models. This presents and interest since allows 
a easier treatment of boundary conditions. 

To simulate the stochastic process, a numerical scheme is proposed and dis- 
cussed into details. This numerical scheme has been developed with some fun- 
damental constraints in mind: 

(i) The numerical scheme must he explicit, uncoditionnaly stable, of order 2 
in time and the number of calls to particle localization has to be minimum. 

(ii) The numerical scheme must be consistent with the analytical solutions of 
the system when the coefficients are constant. 

(iii) The numerical scheme must be consistent with all limit systems. 

These constraints are not mathematical strangeness but they come out from the 
intrinsic structure of the stochastic system and are important for practical con- 
cerns. Indeed, they assure consistency, accuracy and efficiency to the numerical 
method, even for equations which present a multi-scale character. In particular, 
in bounded flows velocity time-scale goes to zero with approaching to the wall 
and the stochastic system (|15p becomes stiff. An algorithm which is not consis- 
tent with this asymptotic limit and it is not uncoditionnaly stable would require 
a related small integration time-step, making numerical simulations practically 
impossible. The consistency has been demonstrated analytically. 

The numerical method obtained is validated in a Hybrid Eulerian/Lagrangian 
configuration. Usually a new numerical scheme is validated with a study of dif- 
ferent errors arising from the numerical scheme, in a stand-alone configuration. 
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Nevertheless, the problem has been analyzed exhaustively in the case of free- 
shear fluid flows with another weak second-order scheme |4] and in the case of 
two-phase flows for an analogous scheme [13]. In particular, the last scheme is 
retrieved rigorously by the present scheme in the case of high Reynolds-number 
flows (when viscosity is put to zero). Thus, we have preferred to concentrate our- 
selves on the major point of the consistency of the hybrid Eulerian/Lagrangian 
method, which is a fundamental issue for this kind of approach [TJ. The global 
method is found to be perfectly consistent at the level of first two velocity 
moments. Moreover, analytical results for Reynolds stress are recovered in the 
asymptotic limits, which supports also the numerical validation of the numerical 
scheme. 

In this numerical validation, it has been pointed out two other main points: 
(i) Eulerian/Lagrangian methods which are not consistent with the same tur- 
bulence model like DNS/SLM or V2F/SLM suffer from an important bias error 
which make result flawed. Therefore such configurations often proposed in two- 
phase flows should reconsidered critically, (ii) It has been shown that in a hybrid 
framework the viscous and non-viscous models give basically the same results, 
provided the correct boundary conditions are imposed. 

Finally, the main point to underline, about the development of the numeri- 
cal scheme, is that the methodology presented goes beyond the borders of the 
present model and of present calculations. The main objective was to propose 
a safe guideline to follow, when one has to deal with numerical simulations of 
stochastic models. Moreover, the methodology is not restricted to fluid mechan- 
ical applications, but it can be applied without modifications to whatever model 
which has a form of linear SDEs, for example in polymeric fluids [1]. 



A Properties of the physical model 

The mean momentum equation is simply obtained by applying the averaging 
operator to the particle velocity equation ([6]) 

Using the relation between the instantaneous substantial derivative and its 
Eulerian counterpart, d ■ /dt = d ■ /dt + Ujd ■ /dxj, we obtain 

d{U.) ^ ,jj,dm I d{u,u,) _ IdjP) ^ ^ d^m ^ .^2) 
dt ^ dxj dxj p dxi dxl 

Thus, the exact mean Navier-Stokes equation is satisfied. This should not be too 
surprising since convection is treated without approximation by the Lagrangian 
point of view and since the mean pressure-gradient, which represents the mean 
value of the acceleration of a fluid particle, is properly taken into account in 
Eq. ([6]). For the second order, one has to write the instantaneous equations for 
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the fluctuating velocity components along a particle trajectory. This is done by 
writing Ui = Ui — (Ui) and consequently 

dt dt \ dt ^ ■" dxj J ^ dxj ^ ' 

We now write the equation in an incremental form to properly handle the 
stochastic terms, which using the mean Navier-Stokes equation is 

du, = dt - uk^^ dt + G^kUk dt + A^kUk dt + ^/Ca{e)dW^. (A.4) 

axk oxk 

The first two terms on the rhs are exact and are independent of the form of 
the stochastic model. The different SDEs are defined in the Ito sense and the 
derivatives of the products UiUj are obtained from Ito's formula 

d{uiUj) = Uiduj + Ujdui + Co{e) dtSij. (A. 5) 

The mean second-order equations are then 



(Uk) + = - ("'^fc)^^ - ("J"^)- 



dt dxk dxk dxk dxk 

+ Gik{ujUk) + Gjk{uiUk) (A. 6) 

^ d{uiUj) 

4' 



B Calculus of the stochastic integrals 

Here, it is explained how the stochastic integrals, appearing in the analytical 
solutions of the equation system with constant coefficients, can be re-arranged 
(by stochastic integration by parts), to yield the covariance matrix. 



B.l Integration by parts 

Let X{t) and Y{t) be two diffusion processes. One can show that (see for 
example Klebaner |12|). in the Ito sense, 

X{t)Y{t)-X{to)Yito)= f Xis)dYis)+ j Y{s)dX{s) + [X,Y]{t), 

J to J to 

where [X, y](t) is the quadratic covariation of X{t) and Y(t) on [tQ,t]. In the 
case where one of the processes is deterministic, [X, Y]{t) = 0. In the frame of 
our study, where the integrated variable is always a deterministic function, one 
can therefore apply integration by parts as in classical differential calculus. 

In fact, in the analytical solutions of the equation system with constant 
coefficients, one encounters multiple stochastic integrals of the type 

1= [ exp{-s/a)( [ exp{s'/b)dW{s')] ds, (B.l) 

Jto \Jto J 
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where (a, 6) G M+^. By setting 

F(s) = [ exp{s'/b)dWis') =^ dF{s) =exp{s/b)dW{s), 

J to 

dG{s) — exp(— s/a) G{s) — —a exp(— s/a) ds, 
and applying integration by parts, one obtains 



/ = — a exp(— </a) / eyip{s/b) dW{s) + a / exp(— s/a) exp(s/6) dVF(s). 

(B.2) 

Therefore, by stochastic integration by parts, the multiple integral given by Eq. 
(jB.ip can be written as the sum of two simple stochastic integrals, Eq. (|B.2p . 



B.2 Derivation of the covariance matrix 

By using the results of the previous subsection and the main properties of the 
Ito integral, that is, linearity, the zero mean property, 

t 

X{s)dW{s)) = 0, 

to 

and the isometry property, 

X{s)dW{s) [ ' Y{s)dW(s)) = [ \x{s)Y{s)) ds, 

Jt2 Jt2 

with to < t2 < ti < t-^. From the zero mean property, it follows that the first 
order moments are equal to zero. For the second order moments (covariance 
matrix), the previous properties give the following equality 

{(j29rn{t)j'u{s)dW{s)^ ) = 

E^mW / f^{s)ds + 2Y,9m{t)gk{t) f f,n{s)fk{s)ds. 

m "'*o "'to 

(B.3) 

where gim{t) and fim(t) are deterministic functions of time. 



C Simulation of a Gaussian vector 

Let X = [Xi, . . . , Xd) be a Gaussian vector defined by a zero mean and a 
covariance matrix Cij = (XiXj). For all positive symmetric matrix (such as 
Cij), there exists a (low or high) triangular matrix Pij which verifies 

d 

fc=i 
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P is given by the Choleski algorithm (here for the low triangular matrix) 



Cil 



1 s^ii^d 




) 



1 <i ^ d 



1 < j < i ^ d 



Pij ==0, i < j d. 



Let G — (Gi, . . . , Gd) be a vector composed of independent Af{0, 1) Gaussian 
random variables, then it can be shown that the vector Y = PG is a Gaussian 
vector of zero mean and whose covariance matrix is C = PP*. Therefore, X 
and Y are identical, that is, 



Eq. (jG.l|) shows how the stochastic integrals, obtained in the analytical solutions 
of the system with constant coefficients, can be simulated. 
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Table 4: Weak semnH order snheme: snh2c 



Prediction step: 

apply the Euler scheme Schl, see Table [3J 



Correction step: 



n+l 



1 



1 



= - {/« exp(-At/7^") + - eM-^t/t2Tr') 



• A(At, rn [Tpcn + B{At, t2T^+') [t2Tl'+H2C^+'] + tsTf^' 



The coefBcients A, B, Ac, Be et Cc are defined as 



A{At, X) ^ ~ exp(-At/X) 



1 - cxp(-A</X) 



At/X 



B{At,X) = 1 



1 - exp(-AVX) 



At/X 



Ac{X,Y) = - exp(-At/X) + 



X + Y 
At 



[1 - exp(-Ai/X)] 



B,iX,Y) = 1 - rl±l[l -exp(-At/X)] + |-c,(x,r), 



CeiX,Y) 



Y 



Y - X 



[exp(-Ai/i") - exp(-Ai/>^)]. 



1 + 1^) Cc{X,Y), 



The stochastic integrals t27i"^^ ^^"^ i2r"^^ are simulated as follows 



7, 



,n+l 



KA^(2At,t2T;"+0 + ar'S7(2At,t2Tf+^)]/i,,(*2T, 



i={2,3} 



n+l 



i=l 



7r * - KA(2At,t2rr') + ar'B^{2At,t2Tr )]hA{t2T, 
+ t2A"+it2r+iKA(2Ai,t2r2"+i) + a^'+iB^(2Ai,t2T2"+i)]/i,2(t2r2"+i: 
- i2^"+'t2r+V^'v4^(2At,t2Ti"+i) + a^+^B^{2At,t2T^+^)]l2,2{t2Tp+^) 



A-y{At,X) = 
B-y{At,X) = 



-cxp{-2At/X) 



1 - exp(-2At/X) 



2 At/X 



1 



1 - exp(-2At/X) 



1 - 



1 - cxp(-2Ai/^) 
2AV^ 



1 - cxp(-2Af/X) 
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